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Abstract 

We  report  on  the  isogeometric  residual-based  variational  multiscale  (VMS)  Large  Eddy 
Simulation  of  a  fully  developed  turbulent  flow  over  a  wavy  wall.  To  assess  the 
predictive  capability  of  the  VMS  modeling  framework,  we  compare  its  predictions 
against  the  results  from  direct  numerical  simulation  (DNS)  and,  when  available, 

against  experimental  measurements.  We  use  C1  quadratic  B-spline  basis  functions  to 
represent  the  smooth  geometry  of  the  sinusoidal  lower  wall  and  the  solution  variables. 
The  Reynolds  number  of  the  flow  considered  is  6,760  based  on  the  bulk  velocity  and 
average  channel  height.  The  ratio  of  amplitude  to  wavelength  ( a/2 )  of  the  sinusoidal 
wavy  surface  is  set  to  0.05.  The  computational  domain  is  22x1.052x2  in  the 
streamwise,  wall-normal  and  spanwise  directions,  respectively.  Mean  averaged 

quantities,  including  velocity  and  pressure  profiles,  and  the  separation/reattachment 
points  in  the  recirculation  region,  are  compared  with  DNS  and  experimental  data.  The 
turbulent  kinetic  energy  and  Reynolds  stress  are  in  good  agreement  with  benchmark 
data.  Coherent  structures  over  the  wavy  wall  are  observed  in  isosurfaces  of  the  Q- 
criterion  and  show  similar  features  to  those  previously  reported  in  the  literature. 
Comparable  accuracy  to  DNS  solutions  is  obtained  with  at  least  one  order  of 

magnitude  fewer  degrees  of  freedom. 

Keywords:  Variational  multiscale  modeling,  Large  Eddy  Simulation,  Wavy  wall, 
Isogeometric  analysis,  B-spline  finite  elements 

1.  Introduction 

Fully  developed  turbulent  flows  over  wavy  surfaces  are  relevant  for  many  engineering 
applications,  such  as,  pipe  flow  in  a  heat  exchanger  where  wall  waviness  enhances 
heat  transfer  and  gas-liquid  contractors  in  the  chemical  industry.  The  wavy  wall 
induces  the  alternation  of  favorable  and  adverse  pressure  gradients  and  a  large  ratio  of 
amplitude  to  wavelength  (a/2)  generates  a  complex  recirculating  flow,  which  enhances 
turbulent  mixing. 

Many  experimental  and  numerical  studies  have  been  conducted  to  understand  the 
complex  physics  of  fully  developed  turbulent  flows  over  wavy  walls.  Hudson  et  al. 
[1]  measured  the  spatial  and  temporal  variation  of  the  velocity  components  in  a  flow 
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with  wavy  walls  using  Laser  Doppler  Velocimetry  (LDV)  and  reported  mean  velocity 
fields  and  fluctuating  quantities  such  as  turbulent  stresses  and  kinetic-energy  production. 
In  their  experiment,  the  Reynolds  number  based  on  bulk  velocity  and  channel  height 
is  6,760  and  the  amplitude  to  wavelength  ratio  is  equal  to  0.05.  These  measurements 
are  used  to  benchmark  our  simulations.  The  experiments  showed  that  the  shear  layer 
in  the  vicinity  of  the  separation  region  is  the  locus  of  large  turbulent  kinetic  energy 
production  and  Reynolds  stresses. 

Gunther  and  Rohr  [2]  studied  streamwise  structures  and  the  three-dimensional  flow 
field  using  proper  orthogonal  decomposition  (POD)  in  their  particle  image  velocimetry 
(PIV)  experiment.  They  identified  the  characteristic  length-scale  of  the  longitudinal 
structures  as  1.52  in  both  laminar  and  turbulent  wavy  wall  flows  and  POD  results 
revealed  smaller  structures  at  the  maximum  Reynolds  stress  location. 

Maass  and  Schumann  [3]  performed  direct  numerical  simulation  (DNS)  of  turbulent 
wavy  wall  flows  using  the  finite  difference  method  and  showed  the  effective  friction 
velocity  is  about  50%  larger  at  the  wavy  surface  than  at  the  flat  surface  because  of 
the  additional  pressure  drag. 

Cherukat  et  al.  [4]  conducted  DNS  studies  of  the  turbulent  flow  over  a  sinusoidal 
wavy  surface  using  spectral  elements  of  7th  order.  Their  results  were  validated  with 
the  experimental  data  of  Hudson  et  al.  [1]  and  are  considered  a  reliable  benchmark 
for  numerical  simulation  of  wavy  walls.  They  report  instantaneous  flow  fields  and 
turbulent  statistics  and  used  them  to  find  velocity  bursts  in  the  separated-flow  region 
that  extend  over  large  distances  away  from  the  wall. 

Recently,  Yoon  et  al.  [5]  investigated  the  effect  of  wave  amplitude  on  turbulent  flow 
in  a  wavy  channel  using  DNS  with  the  immersed  boundary  method.  They  observed 
that  the  pressure  drag  coefficient  increases  as  the  ratio  a/1  increases,  whereas  the 
friction  drag  coefficient  is  maximum  at  the  specified  amplitude  to  wavelength  ratio, 
a/ X=0.3. 

In  [6],  Park  et  al.  compared  several  Reynolds-Averaged  Navier  Stokes  (RANS)  models, 
such  as  standard  k-s  and  linear/nonlinear  k-c-fM  and  showed  that  most  RANS  models 
tested  have  severe  limitations  in  predicting  flow  separation  appropriately  with  the 
exception  of  the  nonlinear  k-s-fM  model.  The  nonlinear  k-E-fM  model  shows  good 
agreement  with  DNS  results  in  mean  flow  field,  streamwise  fluctuations,  and  Reynolds 
stresses,  however,  it  exhibits  discrepancies  in  other  turbulent  quantities  such  as  wall- 
normal  fluctuations. 

Most  of  the  previous  numerical  simulations  of  the  fully  developed  turbulent  flow  over 
a  wavy  wall  have  been  done  using  DNS  to  capture  the  separated  flow  field  and  to 
collect  turbulent  statistics.  However,  very  fine  resolutions  are  required  to  resolve  the 
Kolmogorov  length  scale  of  the  boundary  layer,  which  implies  a  large  computational 
burden.  To  overcome  these  difficulties,  we  believe  that  Large  Eddy  Simulation  (LES) 
is  an  adequate  tool,  since  separation  and  turbulence  fluctuations  can  be  captured  well 
while  utilizing  significantly  less  computational  resources. 

The  recently  proposed  Variational  Multiscale  Method  [7,8]  has  proven  to  be  an 
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effective  approach  to  LES  and  yields  new  insights  into  classical  LES  challenges  such 
as  scale  separation  and  closure  modeling.  The  proposed  residual-based  variational 
multiscale  (VMS)  framework  uses  analytical  tools  to  approximate  the  effect  of  the  fine 
scales  in  the  coarse-scale  governing  equations,  such  as,  perturbation  series  arguments 
and  the  fine-scale  Green’s  function;  see  e.g.  [7].  In  this  framework,  the  simplest 
approximation  to  the  fine  scales  is  given  in  the  form  of  a  scaled  residual  of  the 
coarse  scales.  This  VMS  methodology  was  validated  in  bypass  transition  [9],  forced 
isotropic  turbulence  and  fully-developed  channel  flows  [7,10],  on  free- surface  [11], 
Taylor-Couette  flows  [12],  and,  in  conjunction  with  the  weak  imposition  of  Dirichlet 
boundary  conditions,  in  high  Reynolds  number  channel  and  asymmetric  diffuser  flows 

[13]. 

In  this  study,  we  further  validate  isogeometric  VMS  methodology  in  a  turbulent  flow 
with  separation,  that  is,  in  fully-developed  turbulent  flow  over  a  wavy  wall.  We  report 
mean  fields  and  turbulent  statistics  of  the  flow  such  as  Reynolds  stress,  kinetic  energy 
and  production.  Additionally,  instantaneous  velocity  fields  are  studied  and  their 
coherent  structures  described. 

The  organization  of  the  paper  is  as  follows.  In  Section  2,  following  [7],  we  briefly 
summarize  the  residual-based  variational  multiscale  formulation  for  the  incompressible 
Navier-Stokes  equations  and  introduce  the  numerical  procedure  used  to  discretize  the 
resulting  non-linear  system  of  equations.  In  Section  3,  we  present  the  problem  setup 
describing  the  computational  domain,  boundary  conditions  and  meshes  used.  In  Section 
4,  we  present  the  main  results  for  the  wavy  wall  simulations  and  compare  these  with 
DNS  results  and  experimental  data.  In  section  5  we  draw  conclusions. 

2.  Residual  based  Variational  Multiscale  Method 

The  incompressible  Navier-Stokes  equations  are: 


- \-V -(u®u)  +  Vp  =  vAu  + f  in  f2  (1) 

dt 

V  h  =  0  in  Q  (2) 

w(b+ )—  w(0  )  in  (3) 

»  =  0  on  T  (4) 


where  u  is  the  flow  velocity,  p  is  the  pressure,  /: f2— >R3  is  a  given  body  force,  v  is 
the  kinematic  viscosity,  and  the  symbol,  0  denotes  the  tensor  product.  Equations  (1) 
and  (2)  represent  the  balance  of  linear  momentum  and  mass,  respectively,  while 
equations  (3)  and  (4)  describe  the  initial  and  boundary  conditions,  where  Uo  is  the 
given  initial  velocity  distribution.  Let  V  denote  both  the  trial  solution  and  weighting 
function  spaces,  which  are  assumed  to  be  identical.  The  variational  formulation  can  be 
expressed  as  follows: 

Find  U  ef  such  that  VW  E  V: 


B(W,  U)  =  Bx  (W,  U)  +  B2  (W,  U,  U)  =  L(W ) 
with 
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(5) 


BX(W,U)  = 


+  (q,  V  •  u)n  -  (V  •  w,  p)n  +  (Vsw,2vVsu)n 


(6) 


f  du  ' 
w,— 

V  ^  i 

B2(W,U,V)  =  -{Vw,u®v)a 

L{W)  =  {w,f)n 


(7) 

(8) 


where  1/={m,  p}  and  W={>v,<7}  belong  to  V.  B i(v)  is  a  bilinear  form  and  52(\v)  is  a 
trilinear  form. 


The  solution  and  weighting  function  spaces  are  decomposed  into  coarse-scale  and  fine- 
scale  subspaces  as  follows: 

V  =  V  ©  V  (9) 

where  V  is  the  finite-dimensional  space  associated  with  the  coarse  scales  of  the 
problem  and  V'  is  the  complement  of  V  in  V  and  is  associated  with  fine-scale 

variables.  We  indentify  V  with  the  numerical  solution  of  the  resulting  discrete 
variational  system.  Using  the  linearity  of  [5]  with  respect  to  the  weighting  functions, 
we  divide  the  variational  system  into  a  coupled  coarse-scale  and  fine-scale  equation 
system: 


B§yJj  +  U')=L§y\  WeU 

(10) 

B^V',U  +  U'^)=  l(fVf),  VW’eV' 

(11) 

where  UeV.U’eV 

B§yjj  +  u')=Bx  §y  bx§v  ,u'\  b2§v  jj  jiy  b2§v  jj  ,u) 

+  B2  ({V,  U’.uy  B2  (l V ,  U',  U') 

B$V',U  +  U'y  B \  uyBl  (fVf,  U')+  B2  ([¥’,  U,uyBl  ($¥’,  U,  U’ ) 

(12) 

(13) 

+  B2  U’.uy  B2 (fVr, U\ (/') 

Using  integration  by  parts,  we  rearrange  equation  (13)  as: 

Bv  (W,  U’)+  B,  (W,  V,  £/')=  (w,  Res(t7)^ 

(14) 

where 

B-  (W,  U ')  =  Bx  (VT,  U')+  B2  fy”,  U',  U yB2  ty',  U,  U’) 

(15) 

(w,  Res^)  =  L(W')~  B,  (f,  v\ B,  ^V’.U.u) 

(16) 

Res (U)  represents  the  residual  of  coarse-scales  and  (y)  is  the  duality  pairing 
between  V  and  V'‘ .  Thus,  we  think  of  U'  as  the  solution  of  the  non-linear  problem 
described  by  (14).  It  is  abstractly  expressed  as  a  function  of  U  and  Res (U)  as 
follows: 


U'  =  F'  ^7,  Res^/J 
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(17) 


Following  [7],  once  we  introduce  a  disctretization  and  indentify  the  coarse  scales  with 
the  resolved  scales,  Uh,  we  approximate  the  fine-scale  Green’s  operator,  which  arises 
in  the  perturbation  series  solution  of  (15),  with  an  element-wise  stabilization  operator 
r  and  residual  of  the  coarse  scales  Rqs(Uh )  as  follows: 

U'~  -rRes(Uh)  (18) 

These  approximate  fine  scales  are  expressed  in  terms  of  the  resolved  velocity  and 
pressure  fields  as 


’Ph} 

-rcrc(«")  J 


(19) 


where  fy/  and  rc  are  the  coarse-scale  residuals  of  the  momentum  and  continuity 
equations,  respectively.  The  scaling  factors,  also  known  as  intrinsic  times,  t M  and  tc  , 

are  element-wise  stabilization  parameters  for  the  momentum  and  continuity  equations. 
The  residuals  and  scaling  factors  are  defined  below: 


rM(u",ph)  =  ^  +  uh  •  V«"  +  V/  -vA uh-f 


dt 


rc(uh )  =  V  •  uh 


t  M 


4  -  +  uh  Guh  +Ctv2G:G 


-1/2 


At 2 


■  S)"1 

G  =  f  MtMr 

*  k  a*(  a*, 


§t 


(20) 


where  x.  is  a  physical-space  coordinate  and  is  a  parameter- space  coordinate  (i.e., 

reference  element  coordinate).  Ci  is  set  to  36  through  derivation  of  an  inverse  estimate 
[7].  Finally,  the  formulation  of  the  residual-based  variational  multiscale  method  reads 
as  follows: 


Find  Uh  such  that  \/Wh: 

BM\Wh,Uh)-LM\Wh)  =  0  (21) 

Bm y  bg  tyh  v h  ) 

+  («*  ■Vi»‘  +  V,r„r„  («*,/>*  1  +  (v-^‘,Tcrc («")),  (22) 

+  («*,/>*))-  ,TMrM(fih  ,ph)ai:uru(fih  ,p‘‘\ 
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(23) 


LM\Wh)  =  (w\f), 


Bu(W\Uh)  = 


„  3 «*  ' 

>V 


3/ 


-<V» ^«‘®«*M(v-»’V)0+(VV,2l/VV), 


(24) 


+(*\W)fl 

where  the  superscripts  G  and  MS'  stand  for  Galerkin  and  variational  multiscale, 
respectively. 


We  use  the  generalized-a  method  as  presented  by  Jansen  et  al.  [14]  to  integrate  the 
governing  equations  in  time.  We  set  p=0.5,  where  p,„  is  the  spectral  radius  of  the 
amplification  matrix  as  f— ><= «.  To  solve  the  nonlinear  system  of  equations,  we  employ 
Newton's  method  with  a  two-stage  predictor-multicorrector  algorithm.  See  [7]  for 
further  details. 


3.  Numerical  simulation 


The  Reynolds  number  is  6,760  based  on  the  bulk  velocity  and  average  channel  height. 
The  geometry  of  the  wavy  wall  is  described  by  the  formula  acos(2nx/l),  where  a  is 
the  amplitude  of  the  wave  and  2  is  the  wavelength.  Herein,  the  ratio  (a/2)  of  the 
amplitude  to  the  wavelength  is  set  to  0.05,  which  is  the  same  as  that  in  references 
[3-5], 


The  computational  domain  ( LxxLyxL:)  shown  in  Figure  1  has  dimensions  of 
22x1.052x2  in  the  streamwise  (v),  vertical  (y)  and  spanwise  (z)  directions.  This  is 
identical  to  the  one  used  by  Yoon  et  al.  [5].  We  employ  no-slip  Dirichlet  boundary 
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conditions  at  the  lower  wavy  wall  and  at  the  upper  flat  wall  (y=l .()/),  and  we  apply 
periodic  boundary  conditions  in  both  the  streamwise  and  spanwise  directions.  The 
driving  force  is  the  streamwise  pressure  gradient  which  is  adjusted  to  maintain  a 
constant  mass  flow  rate. 

Isogeometric  analysis  [15]  with  B-spline  basis  functions  is  employed  in  the  present 
work.  The  control  points  for  the  B-spline  curve  defining  the  geometry  are  obtained  by 
projecting  the  known  sinusoidal  geometry  onto  the  coarsest  possible  quadratic  B-spline 
space  that  can  represent  the  geometry  [16,  17].  The  open,  non-uniform  knot  vector 
is: 


S  =  {0, 0,0, 1,2, 3, 4, 5, 6, 7, 8,8, 8} 


and  the  control  points  are  reported  in  Table  1. 


Table  1.  Control  points  for  the  B-spline  curve  defining  the  bottom  wavy  wall. 


X 

y 

5, 

0.0000E+00 

5.0000E-02 

b2 

3.1250E-02 

4.9998E-02 

9.3750E-02 

4.2370E-02 

b4 

1.5625E-01 

2.8299E-02 

b5 

2.1875E-01 

9.9351E-03 

b6 

2.8125E-01 

-9.9350E-03 

Bn 

3.4375E-01 

-2.8299E-02 

B 8 

4.0625E-01 

-4.2370E-02 

B 9 

4.6875E-01 

-4.9998E-02 

B\q 

5.0000E-01 

-5.0000E-02 

4.  Results:  Simulation  of  wavy  wall  channel  flows 

We  analyzed  the  sensitivity  of  the  simulation  results  to  the  grid  resolution  using  the 
following  three  cases:  128x64x64  (Fine  mesh),  64x32x32  (Medium  mesh)  and 
32x16x16  (Coarse  mesh).  The  resolutions  correspond  to  the  streamwise,  vertical  and 
spanwise  directions,  respectively.  For  all  meshes,  ^-continuous  quadratic  B-splines 
basis  functions  are  used.  The  geometry  definition  and  solution  variables  use  identical 
parameterization  (isoparametric  disctretization).  The  number  of  mesh  elements  for  the 
three  cases  is  summarized  and  compared  with  meshes  used  in  other  DNS 
computations  in  Table  2.  An  equivalent  number  of  mesh  points  per  unit  streamwise 
and  unit  spanwise  wave  length  are  reported  in  Table  2  to  allow  a  fair  comparison  of 
the  different  simulations,  since  different  computational  domain  sizes  were  used  in  each 
simulation.  As  mentioned  previously,  the  spectral  element  method  of  Cherukat  el  al. 
[4]  employs  7th  order  Lagrange  polynomials.  The  total  number  of  mesh  points  is 
64x148x32,  which  is  greater  than  3.0xl05.  In  the  other  two  DNS  cases,  based  on 
finite  difference  and  finite  volume  methods,  3.0xl05  and  1.8xl06  mesh  points  are  used. 
The  number  of  degrees  of  freedoms  in  the  medium  mesh  case  of  the  present  work  is 
3.3xl04  (=323),  which  is  at  least  an  order  of  magnitude  fewer  than  those  used  in  the 
DNS  calculations. 

Table  2.  Mesh  elements  in  the  present  work  and  in  the  DNS  calculations  of  others. 
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Lx 

Ly 

Lz 

Nx 

Ny 

Nz 

NA 

NA 

Coarse  (present) 

32 

16 

16 

16 

16 

Medium  (present) 

22 

1.052 

2 

64 

32 

32 

32 

32 

Fine  (present) 

128 

64 

64 

64 

64 

Cherukat  et  al.  [4]* 

42 

1.052 

22 

36 

21 

64 

9 

32 

Maass  and  Schumann  [3] 

42 

1.052 

22 

256 

128 

96 

64 

48 

Yoon  et  al.  [5] 

22 

1.052 

2 

250 

150 

100 

125 

100 

Spectral  method  with  7th  ore 

er  Lagrange  po 

ynomials. 

Figure  2.  Grid  convergence  test,  top:  pressure  coefficient  distribution  at  the  wall,  bottom:  urms,  vrms  at 

two  positions  (a)  x/2=0.0  (crest),  (b)  x/A=0.5  (trough). 


Results  for  the  three  meshes  are  presented  in  Figure  2.  The  mean  pressure  coefficient 
at  the  wavy  wall  shows  that  the  coarse-mesh  result  has  small  discrepancies  compared 
with  the  two  finer  cases.  The  root-mean- square  (rms)  velocity  fluctuations  are 
compared  in  Figure  2  at  two  positions:  the  crest  (x//=0.0)  and  the  trough  (x//=0.5). 
Velocity  fluctuations  of  the  medium  and  fine  meshes  are  consistent  with  each  other, 
while  the  coarse  mesh  fails  to  properly  resolve  the  velocity  fluctuations.  This  lack  of 
resolution  is  particularly  noticeable  in  the  region  of  rapid  variation  near  the  peak. 
Based  on  these  results,  the  medium  mesh  (64x32x32)  is  considered  to  be  adequate. 
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To  compute  the  mean  flow  quantities,  we  averaged  the  instantaneous  velocity  and 
pressure  fields  in  the  spanwise  direction  and  in  time  over  a  period  of  100 H/Ub,  where 
H  is  the  channel  height  and  Ui,  is  the  mean  bulk  velocity. 


Figure  3  shows  mean  velocity  vectors  and  streamlines  for  the  medium  mesh  case.  It 
can  be  seen  that  the  turbulent  boundary  layer  near  the  wall  is  well  resolved. 


Figure  4.  Shear  stress  at  the  wall  (medium  mesh  case). 

To  determine  the  separation  and  reattachment  points,  the  wall  shear  stress  is  calculated 
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and  shown  in  Figure  4.  The  separation  point  is  seen  to  be  predicted  at  x/}.-0. 1 5  and 
the  reattachment  point  at  x//=0.58.  The  DNS  results  of  Cherukat  et  al.  [4]  give  the 
corresponding  points  at  0.14  and  0.59,  and  the  experiment  of  Hudson  et  al.  [1] 
determined  these  points  as  0.22  and  0.58.  As  argued  by  Cherukat  et  al.  [4],  this 
difference  in  the  measured  separation  point  between  numerical  simulation  and 
experiment  may  be  due  to  the  difficulty  of  measuring  separation  when  there  are  large 
oscillations  in  the  velocity  while,  at  the  same  time,  the  mean  flow  velocity  is  small. 
The  prediction  of  these  points  as  0.15  and  0.59  in  the  DNS  of  Maass  and  Schumann 

[3]  is  also  consistent  with  the  present  results. 

The  shear  stress  increases  rapidly  after  the  minimum  value  near  x/2=0.40,  and  has 
maximum  value  near  x/2=0.92,  which  is  in  agreement  with  reference  data  reported  in 

[4] .  The  shear  stress  at  the  upper  flat  plate  wall  is  computed  as  0.0043  in  the  present 
work,  compared  with  the  value  0.0045  of  [4]. 


Figure  5.  Mean  streamwise  velocity  profiles,  (a)  x/X= 0.1,  (b)  x/2=0.3,  (c)  x/l=0.5,  (d)  x//=().l  (medium 

mesh  case). 
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Figure  6.  Mean  vertical  velocity  profiles,  (a)  x/X= 0.1,  (b)  x//=0.3,  (c)  x//=0.5,  (d)  x//=().l 

(medium  mesh  case). 


Figures  5  and  6  show  the  streamwise  and  vertical  mean  velocity  profiles  at  different 
streamwise  locations  (xA=0.l,  0.3,  0.5,  0.7)  and  these  are  compared  with  the  DNS 
results  [3].  Point  a  (x/)=().  1 )  is  located  near  the  crest  of  the  wavy  wall  and  before  the 
separation  point  of  the  flow.  Points  b  and  c  (x/A=0. 1 ,  0.5,  respectively)  are  in  the 
recirculation  region.  Point  d  (v//.=0.7)  is  located  after  the  reattachment  point  of  the 
flow.  The  streamwise  and  vertical  velocity  profiles  agree  well  with  the  DNS  results 
reported  in  [3]. 


11 


Figure  7.  Mean  pressure  coefficient  distribution,  line:  present  simulation,  symbol:  DNS  of  Cherukat  et 

al.  [4]  (medium  mesh  case). 


The  mean  pressure  coefficient  along  the  wavy  wall  is  replotted  in  Figure  7  along  with 
the  locations  of  separation  and  reattachment  points.  The  separation  bubbles  are  seen  to 
appear  within  the  adverse  pressure  gradient  regions  as  is  typical.  The  peak  pressure 
point  occurs  near  x/X=0.69.  From  that  point  onward  a  favorable  pressure  gradient 
dominates  until  the  next  crest  point  at  x/X- 1.0. 


Figure  8.  Reynolds  stress  contours,  top:  present  simulation,  bottom:  DNS  of  Cherukat  et  al.  [4] 

(medium  mesh  case). 
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Figure  8  shows  the  contours  of  the  Reynolds  stress  term,  -u'v'.  The  upper  plot 
displays  the  present  results  and  the  lower  one  reproduces  those  of  [4].  As  shown  in 
Figure  8  (top),  the  position  of  maximum  Reynolds  stress  is  predicted  near  the  upper 
region  of  the  trough  (x//.=0.45)  and  negative  values  are  seen  on  the  windward  side 
near  the  wall  between  x/2=0.53  and  0.99.  The  DNS  results  of  Cherukat  et  al.  [4]  give 
the  maximum  value  at  v//=0.40  and  the  negative  region  is  from  x//=0.55  to  x/X=  1.0. 
Hudson  [1]  and  others  [4,  5]  argue  that  the  negative  values  of  Reynolds  stress  occur 
due  to  the  calculation  of  this  term  in  a  Cartesian  frame  rather  than  in  a  boundary 
fitted  frame. 

Gunther  and  Rohr  [2]  described  these  wavy  flows  in  terms  of  three  regions;  (I) 
separation  zone,  (II)  region  of  local  maximum  Reynolds  stress  (-u'v')  and  (III)  region 
of  local  minimum  Reynolds  stress  (- u'v *).  The  contours  of  maximum  Reynolds  stress 
after  the  reattachment  point  increase  in  an  upward  direction,  which  is  represented  by 
the  dashed  lines  in  Figure  8.  The  local  maximum  in  the  Reynolds  stress  contours  is 
approximately  located  at  v=0.08/.  above  the  wall.  This  coincides  with  Gunther  and 
Rohr's  experimental  results  [2]. 


Figure  9.  Reynolds  stress  contours,  left:  present  simulation  results,  right:  DNS  of  Cherukat  et  al.  [4], 
top:  u'u',  middle:  vV,  bottom:  w 'w '  (medium  mesh  case). 


The  contours  of  the  diagonal  components  of  Reynolds  stress  are  plotted  for 

comparison  with  DNS  results  in  Figure  9.  The  locations  of  maximum  streamwise 

fluctuations  are  coincident  with  those  of  maximum  positive  Reynolds  stress  (-u'v').  A 
ridge  appears  over  the  crest  of  the  streamwise  intensity  plot.  However,  the  regions  of 
maximum  vertical  intensity  are  closer  to  the  wall  than  those  of  the  other  two 
variables.  Spanwise  intensity  contours  show  a  different  pattern  from  those  of 

streamwise  and  vertical  intensity.  A  high  intensity  region  between  the  trough  and  the 
crest  indicates  that  the  flow  has  developed  fully  three-dimensional  structures.  These 

terms  contribute  to  the  increase  of  the  turbulent  kinetic  energy  after  the  reattachment 
point,  as  shown  in  Figure  10. 
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The  profiles  of  the  resolved  turbulent  kinetic  energy  at  different  spanwise  locations  are 
shown  in  Figure  10  and  are  compared  with  the  DNS  results  of  Maass  and  Schumann 
[3].  The  present  results  show  very  good  agreement  with  the  DNS  data  even  though 
small  sharp  peaks  are  apparent  at  the  maximum  points  of  Figures  10(b)  and  10(c). 
The  contours  of  turbulent  kinetic  energy  (not  shown  here)  show  the  maximum  kinetic 
energy  point  at  xA=0A5  which  is  the  same  as  the  DNS  results  reported  in  [3]. 

□  Maass  and  Schumann  [3] 


Figure  10  Turbulent  kinetic  energy  profiles,  (a)  x//=().  1 ,  (b)  x//=0.3,  (c)  x/A=0.5,  (d)  x/ /=().! 

(medium  mesh  case). 
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Figure  11.  Turbulence  production  contours,  (a)  results  for  the  medium  mesh,  (b)  results  for  the  fine 
mesh,  (c)  experiment  of  Hudson  et  al.[l],  (d)  DNS  of  Cherukat  et  al.  [4]. 


The  contours  of  the  production  of  turbulent  kinetic  energy  {Pk=-u'iu 


are 


shown  in  Figure  11.  The  overall  patterns  displayed  by  the  present  simulations  are  in 
reasonable  agreement  with  the  reference  data.  Regions  where  negative  production 
occurs  are  predicted  near  the  bottom  wall,  above  the  high  intensity  region,  between 
the  reattachment  point  and  the  next  crest.  These  negative  production  regions  are 
associated  with  energy  backscatter,  that  is,  transfer  of  energy  from  the  smaller  flow 
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scales  to  the  larger  ones. 


The  regions  of  high  kinetic-energy  production  show  slightly  different  patterns  in 
Hudson's  experiment  [1]  and  in  Cherukat  et  al.  DNS  [4].  Hudson's  experiment  predicts 
this  region  near  the  crest,  while  numerical  simulations  of  Cherukat  et  al.  show  the 
high  intensity  region  located  after  the  reattachment  point  to  the  next  crest.  Cherukat  et 
al.  posited  this  difference  may  be  caused  by  improper  resolution  of  Hudson's 
experiment  near  the  wall  in  the  recirculation  zone,  which  decreases  the  accuracy  of 
the  calculated  velocity  gradient  used  in  the  computation  of  the  production  terms  from 
the  experimental  data.  The  patterns  of  high  value  in  this  region  in  the  present 
calculations  are  qualitatively  consistent  with  the  DNS  results.  Angelis  et  al.  [18] 
explain  the  importance  of  the  gradient  of  spanwise  velocity  in  wavy  wall  flows  and, 
in  particular,  near  the  wall  to  accurately  predict  the  turbulent  kinetic  energy  production 
in  the  flow.  This  is  evident  when  comparing  medium  and  fine  meshes  in  Figures 
11(a)  and  11(b).  In  the  medium  resolution  mesh  simulation,  the  high  intensity  region 
is  located  near  the  crest  of  the  wavy  wall  (x/X  >  0.85),  while,  in  the  fine  mesh 
simulation,  the  high  intensity  region  is  located  toward  the  reattachment  point  (x/X  > 
0.75)  and  is  spread  more  widely  than  for  the  medium  mesh. 


15.0 
12.0 
9.0 
6.0 
3.0 
0.0 
-3.0 
-6.0 
-9.0 
-12.0 

0  0.2  0.4  ,,  0.6  0.8  1  '15  ° 

x/  X 

Figure  12.  (Top)  Isosurfaces  of  coherent  structures  at  instantaneous  time  (Q-criterion)  of  fine  mesh, 
contoured  by  streamwise  vorticity  (cox),  (Bottom)  Spanwise  vorticity  contours  in  x-y  plane  at 

z/X= 0.5. 
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To  study  the  coherent  structures  of  this  fully  turbulent  flow,  the  second  invariant  of 
the  velocity  gradient  tensor,  Q-criterion  [19],  is  calculated  for  the  fine  mesh  simulation. 
Even  though  the  predictions  of  the  mean  and  turbulent  statistics  for  the  medium  mesh 
are  in  excellent  agreement  with  the  benchmark  experimental  and  simulated  data,  the 
medium  mesh  resolution  is  not  able  to  accurately  resolve  the  turbulent  coherent 
structures.  Figure  12  (Top)  shows  the  instantaneous  vortical  structures  represented  by 
the  isosurfaces  of  the  Q-criterion  [19]  which  are  colored  with  the  values  of 
streamwise  vorticity  (ojx).  Figure  12  shows  that  the  dominant  structures  are  streamwise 
vortices;  some  structures  begin  after  the  trough  and  extend  to  the  next  crest  or  trough, 
while  others  are  generated  in  the  region  between  the  reattachment  point  and  the  next 
crest.  These  structures  appear  to  scale  like  the  wavelength  of  the  sinusoidal  wavy  wall. 
Hairpin  vortices  are  found  near  the  reattachment  region  and,  at  this  instant,  one  is 
marked  by  a  circle  in  Figure  12  (Top).  Vortices  aligned  with  the  spanwise  direction 
are  rarely  found  due  to  the  dominant  streamwise  convection  of  vortices.  Figure  12 
(Bottom)  shows  the  contours  of  the  spanwise  vorticity  (coz)  in  the  midsection  (z//=0.5) 
of  the  spanwise  direction.  High  vorticity  is  generated  by  the  shear  layer  interaction 
with  the  recirculating  flow.  This  high  vorticity  moves  toward  the  next  crest  while 
breaking  up.  A  thin  vortex  sheet  is  developed  near  the  wall  after  the  reattachment 
point  and  results  in  the  shear  layer  at  the  next  crest.  These  observations  are  consistent 
with  those  of  other  numerical  simulations  [18,  20]. 

5.  Conclusions 

Fully  developed  turbulent  flows  with  wavy  walls  are  simulated  using  an  isogeometric 
residual-based  variational  multiscale  method.  C1  quadratic  B-spline  basis  functions  are 
used  to  solve  the  incompressible  Navier-Stokes  equations  and  represent  the  geometry. 
To  assess  the  predictive  capabilities  of  the  VMS  modeling  framework,  mean  and 
turbulent  statistics  as  well  as  instantaneous  flow  features  are  compared  with  the  DNS 
results  of  Maass  and  Schumann[3],  Cherukat  et  al.  [4]  and  experimental  data  of 
Hudson  et  al.  [1]. 

Mean  velocity  and  pressure  fields  show  good  agreement  with  DNS  results.  Turbulent 
statistics  at  the  wall,  such  as  wall  shear  stress  and  pressure  coefficient,  are  in  good 
agreement  with  reference  data.  The  separation  point  is  predicted  at  x//.=0. 1 5  and  the 
reattachment  point  at  a'//.=0.58.  in  good  agreement  with  the  estimation  of  the  size  of 
the  recirculation  bubble  reported  in  the  literature  [4]. 

Overall,  the  patterns  of  the  contours  of  the  kinetic-energy  production  are  consistent 
with  the  DNS  results.  Additionally,  energy  backscatter  is  correctly  predicted.  As  the 
mesh  is  refined,  particularly  in  the  spanwise  direction,  the  region  of  higher  intensity 
of  production  spreads  toward  the  downwind  side  of  the  wavy  wall,  confirming  the 
observation  of  Angelis  et  al.  [18]  that  the  gradient  resolution  in  the  spanwise  direction 
has  an  important  effect  on  the  calculation  of  production  of  the  turbulent  kinetic  energy. 

Streamwise  vortical  structures  with  a  length  scales  comparable  to  that  of  the  bottom 
wall  wave  length  are  observed  as  dominant  structures  in  this  flow.  Hairpin  structures 
are  also  found  in  this  flow  regime.  Strong  interaction  between  the  shear  layer  and  the 
recirculation  bubble  appears  to  distort  the  thin  boundary  layer  developing  from  the 
previous  reattachment  point. 
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From  the  engineering  point  of  view,  the  residual-based  variational  multiscale  modeling 
framework  is  consistently  able  to  accurately  predict  the  important  features  of  the  mean 
flow  as  well  as  those  of  the  turbulent  statistics.  This  observation  has  been  verified 
herein  for  separated  turbulent  flow  with  relatively  coarse  meshes.  Most  statistical 
quantities  for  the  medium  mesh  are  in  excellent  agreement  with  DNS  data,  but  with  a 
reduction  of  degrees  of  freedom  by  at  least  an  order  of  magnitude.  As  the  mesh  is 
refined,  the  variational  multiscale  results  converge  to  those  of  DNS  results,  as 
anticipated  from  the  structures  of  the  variational  multiscale  modeling  theory. 
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